#Michael Barnett
#4-24-10
#New approach for NAMCS 2005-7 data to examine influence of several patient and 
#physician factors on ambulatory care quality for several cardiovascular disease measures

library(foreign)
library(survey)
	#Adjusting for sample design in survey package, this is why all the warnings don't crash the package
	options(survey.lonely.psu="remove")
	options(survey.ultimate.cluster=TRUE)
library(Zelig)
library(MatchIt)
library(Amelia)
library(cem)
library(ggplot2)

###########
#Section 1: Preparing Database
###########
#Import NAMCS data
load(file = "namcs05.RData")
load(file = "namcs06.RData")
load(file = "namcs07.RData")


#Choose variables in all 3 datasets which will be used in analysis
#Variables for 2005-7
namcs05[, c("pccprod", "pccsat", "pccqoc", "pccpprof", "measpub")] <- -9
vars <- c("age", "sex", "racer", "paytype", "diag1r", "diag2r", "diag3r", "rfv1", "rfv2", "rfv3", "cebvd", "chf", "diabetes", "htn", "ihd", "totchron", "bpsys", "bpdias", "msa", "patwt", "spec", "phycode", "patcode", "retypoff", "solo", "empstat", "owns", "emedrec", "edemog", "ecpoe", "ectoe", "eresult", "epnotes", "eremind", "prprvtr", "year", "physwt", "cstratm", "cpsum", 
"pccprod", "pccsat", "pccqoc", "pccpprof", "measpub",

	"rx1cat1", "rx1cat2", "rx1cat3", "rx1cat4",
	"rx2cat1", "rx2cat2", "rx2cat3", "rx2cat4", 
	"rx3cat1", "rx3cat2", "rx3cat3", "rx3cat4", 
	"rx4cat1", "rx4cat2", "rx4cat3", "rx4cat4", 
	"rx5cat1", "rx5cat2", "rx5cat3", "rx5cat4", 
	"rx6cat1", "rx6cat2", "rx6cat3", "rx6cat4", 
	"rx7cat1", "rx7cat2", "rx7cat3", "rx7cat4", 
	"rx8cat1", "rx8cat2", "rx8cat3", "rx8cat4") 

#Combine three years of data together
	namcs <- rbind(namcs05[,vars], namcs06[,vars], namcs07[,vars] )
	
#Subset to adult patient visits (>18 yo) to primary care and cardiologists;
	namcs <- namcs[namcs$spec %in% c("FPG", "FSM", "GP", "FP", "IM", "IMG", "CD", "ICE") & namcs$age >=18,]
	
#Create a unique ID for each physician and carry physician weight through to all rows
#Unique ID = "physcodeyear"
#	Note: it seems that some physicians don't have a weight ... will exclude them 
	namcs$physcodeyear <- paste(namcs$phycode, namcs$year, sep="")
	physwt_key <- na.omit(namcs[,c("physcodeyear", "physwt")])
	namcs$physwt <- NULL
	namcs <- merge(namcs, physwt_key, by.x="physcodeyear", by.y="physcodeyear", all.y=TRUE)

#Recode variables into workable categories
	
#Physician variables: EHR use and types, office typse, employment, pratice ownership, 
					# specialty, % private insurance patients
					
	namcs$produc <- 0
	namcs[namcs$pccprod==1, "produc"] <- 1
	namcs[namcs$eremind %in% c(-9,-8), "produc"] <- NA # Missing
	namcs$satis <- 0
	namcs[namcs$pccsat==1, "satis"] <- 1
	namcs[namcs$pccsat %in% c(-9,-8), "satis"] <- NA # Missing
	namcs$qoc <- 0
	namcs[namcs$pccqoc==1, "qoc"] <- 1
	namcs[namcs$pccqoc %in% c(-9,-8), "qoc"] <- NA # Missing
	namcs$prof <- 0
	namcs[namcs$pccpprof==1, "pprof"] <- 1
	namcs[namcs$pccpprof %in% c(-9,-8), "pprof"] <- NA # Missing
	namcs$produc <- 0
	namcs[namcs$measpub==1, "public"] <- 1
	namcs[namcs$measpub %in% c(-9,-8), "public"] <- NA # Missing
			
	#Any EHR use
		namcs$EHR <- 0 #No, Don't know, Turned off, or Not applicable
		namcs[namcs$emedrec %in% c(1,2), "EHR"] <- 1 # Yes, part or complete
		namcs[namcs$emedrec %in% c(0,4) & namcs$year %in% c(2005,2006), "EHR"] <- NA # Missing
		namcs[namcs$emedrec %in% c(-9,-8) & namcs$year==2007, "EHR"] <- NA # Missing
		#E-reminder
			namcs$remind <- 0
			namcs[namcs$eremind==1, "remind"] <- 1
			namcs[namcs$eremind %in% c(3,4) & namcs$year==2005, "remind"] <- NA # Missing
			namcs[namcs$eremind %in% c(3,9) & namcs$year==2006, "remind"] <- NA # Missing
			namcs[namcs$eremind %in% c(-9,-8) & namcs$year==2007, "remind"] <- NA # Missing
		#Electronic demographic information
			namcs$demog <- 0
			namcs[namcs$edemog==1, "demog"] <- 1
			namcs[namcs$edemog %in% c(3,4) & namcs$year==2005, "demog"] <- NA # Missing
			namcs[namcs$edemog %in% c(3,9) & namcs$year==2006, "demog"] <- NA # Missing
			namcs[namcs$edemog %in% c(-9,-8) & namcs$year==2007, "demog"] <- NA # Missing		#Electronic prescribing
			namcs$escrip <- 0
			namcs[namcs$ecpoe==1, "escrip"] <- 1
			namcs[namcs$ecpoe %in% c(3,4) & namcs$year==2005, "escrip"] <- NA # Missing
			namcs[namcs$ecpoe %in% c(3,9) & namcs$year==2006, "escrip"] <- NA # Missing
			namcs[namcs$ecpoe %in% c(-9,-8) & namcs$year==2007, "escrip"] <- NA # Missing		#Electronic test orders
			namcs$etest <- 0
			namcs[namcs$ectoe==1, "etest"] <- 1
			namcs[namcs$ectoe %in% c(3,4) & namcs$year==2005, "etest"] <- NA # Missing
			namcs[namcs$ectoe %in% c(3,9) & namcs$year==2006, "etest"] <- NA # Missing
			namcs[namcs$ectoe %in% c(-9,-8) & namcs$year==2007, "etest"] <- NA # Missing
		#Electronic physician notes
			namcs$enotes <- 0
			namcs[namcs$epnotes==1, "enotes"] <- 1
			namcs[namcs$epnotes %in% c(3,4) & namcs$year==2005, "enotes"] <- NA # Missing
			namcs[namcs$epnotes %in% c(3,9) & namcs$year==2006, "enotes"] <- NA # Missing
			namcs[namcs$epnotes %in% c(-9,-8) & namcs$year==2007, "enotes"] <- NA # Missing
		#Electronic test results
			namcs$result <- 0
			namcs[namcs$eresult==1, "result"] <- 1
			namcs[namcs$eresult %in% c(3,4) & namcs$year==2005, "result"] <- NA # Missing
			namcs[namcs$eresult %in% c(3,9) & namcs$year==2006, "result"] <- NA # Missing
			namcs[namcs$eresult %in% c(-9,-8) & namcs$year==2007, "result"] <- NA # Missing
	#Office type
		namcs$office <- 0 #"Other"
		namcs[namcs$retypoff==1,"office"] <- 1 #Private solo or group
		#namcs[namcs$retypoff==2,"office"] <- 2 #Free-standing clinic
		#namcs[namcs$retypoff==3,"office"] <- 3 #CHC
	#Alternate categorization as binary	
		#namcs$office.binary <- 0 #Solo
		#namcs[namcs$retypoff!=1, "office.binary"] <- 1 #Other
	#Employment								        
		namcs$employ <- 0 #Other
		namcs[namcs$empstat==1,"employ"] <- 1 #Owner
		namcs[namcs$empstat==4,"employ"] <- NA #Missing
	#Practice ownership
		namcs$owner <- 0 #Other
		namcs[namcs$owns==1 & namcs$year %in% c(2001:2007),"owner"] <- 1 #Physician
		#namcs[namcs$owns==2 & namcs$year %in% c(2001:2007),"owner"] <- 2 #HMO or Health Care Corporation
	#Percent private insurance
		namcs$private <- namcs$prprvtr
		namcs[namcs$private %in% c(1,2), "private"] <- 0 #<50%
		namcs[namcs$private %in% c(3,4), "private"] <- 1 #>50%
		namcs[namcs$private %in% c(-9,9), "private"] <- NA #Missing
	#Generalist vs. Specialist
		namcs$gp <- 0 #Cardiologist
		namcs[namcs$spec %in% c("FPG", "FSM", "GP", "FP", "IM", "IMG"), "gp"] <- 1 #Generalist
	#Solo practice: recode 2 to 0, 0 = not solo, 1 = solo practice 
		namcs[namcs$solo==2, "solo"] <- 0
	
#Patient variables	
	#Patient payment
		namcs$payment <- 4 #Other or unknown
		namcs[namcs$paytype %in% c(1,4),"payment"] <- 1 #Private, adequate insurance
		namcs[namcs$paytype==2,"payment"] <- 2 #Medicare
		namcs[namcs$paytype %in% c(3,5,6), "payment"] <- 3 #Medicaid, inadequate insurance
		namcs[namcs$paytype %in% c(0,8), "payment"] <- NA #Missing	#Urban/Rural
		namcs$urban <- 0
		namcs[namcs$msa==1, "urban"] <- 1
	#White
		namcs$white <- 0
		namcs[namcs$racer==1, "white"] <- 1
	#Male
		namcs$male <- 0
		namcs[namcs$sex==2, "male"] <- 1
	#Age, in roughly tertiles
		namcs$age.cat <- 1 # <45 years old
		namcs[namcs$age>=45 & namcs$age <65, "age.cat"] <- 2 # 45-64 years old
		namcs[namcs$age>=65, "age.cat"] <- 3 # >=65 years old
	#Number of chronic conditions: 0, 1, 2, 3, 4+
		namcs$numchron <- 0 # 0 conditions
		namcs[namcs$totchron %in% c(-9, 99), "numchron"] <- NA #missing
		namcs[namcs$totchron==1, "numchron"] <- 1 # 1 condition
		namcs[namcs$totchron==2, "numchron"] <- 2 # 2 conditions
		namcs[namcs$totchron==3, "numchron"] <- 3 # 3 conditions 
		namcs[namcs$totchron>=4, "numchron"] <- 4 # 4+ conditions 
		
####################
#Section 2: Functions for data preparation
####################			
diags <- function (codes, type = "diag"){
	#Modified 4/17/10 for NAMCS 2005-7 data
	#Takes a vector of ICD-9, Reason For Visit or Drug codes
	#and returns the logical vector of any NAMCS visit containing an entry for that code
	#enter "diag" for ICD-9 code or "class" for general drug class
	
	if (type == "diag"){
	indices <- sapply(codes, 
					  function (code) 
							c(grepl(code, namcs$diag1r) | grepl(code, namcs$diag2r) | grepl(code, namcs$diag3r)), 
					  simplify=TRUE)
	}
	if (type == "rfv"){
	indices <- sapply(codes, 
					  function (code) 
							c(grepl(code, namcs$rfv1) | grepl(code, namcs$rfv2) | grepl(code, namcs$rfv3)), 
					  simplify=TRUE)	
					  }	
	if (type == "class"){
				indices <- sapply(codes, 
					  function (code) 
							c(grepl(code, namcs$rx1cat1) | grepl(code, namcs$rx1cat2) | grepl(code, namcs$rx1cat3) | grepl(code, namcs$rx1cat4) |
							  grepl(code, namcs$rx2cat1) | grepl(code, namcs$rx2cat2) | grepl(code, namcs$rx2cat3) | grepl(code, namcs$rx2cat4) |
							  grepl(code, namcs$rx3cat1) | grepl(code, namcs$rx3cat2) | grepl(code, namcs$rx3cat3) | grepl(code, namcs$rx3cat4) |
							  grepl(code, namcs$rx4cat1) | grepl(code, namcs$rx4cat2) | grepl(code, namcs$rx4cat3) | grepl(code, namcs$rx4cat4) |
							  grepl(code, namcs$rx5cat1) | grepl(code, namcs$rx5cat2) | grepl(code, namcs$rx5cat3) | grepl(code, namcs$rx5cat4) |
							  grepl(code, namcs$rx6cat1) | grepl(code, namcs$rx6cat2) | grepl(code, namcs$rx6cat3) | grepl(code, namcs$rx6cat4) |
							  grepl(code, namcs$rx7cat1) | grepl(code, namcs$rx7cat2) | grepl(code, namcs$rx7cat3) | grepl(code, namcs$rx7cat4) |
							  grepl(code, namcs$rx8cat1) | grepl(code, namcs$rx8cat2) | grepl(code, namcs$rx8cat3) | grepl(code, namcs$rx8cat4)),
							  simplify=TRUE)
		}	
	#Add up all rows to collect any visit that satisfies criteria
		if(length(indices) != 1) index <- as.logical(rowSums(indices))
		#If only one code, convert to vector
		if(length(indices) == 1) index <- indices[,1]
	print(paste("Number of visits:", sum(index), sep=" "))
	return(index)
}

phys.level <- function(numerator, denominator){
	#Takes as input names of variables from the NAMCS dataset, numerator and denominator
	#Outputs physician level weighted estimate of % difference from average performance
	#for the input quality measure for merging back into data.frame

		print("Splitting data ... ")
		phys <- split(namcs, namcs$physcodeyear)
		avg <- sum(namcs[,numerator])/sum(namcs[,denominator])
			#Center the binary QM measure on the average physician performance
			#e.g: 70% physicians on average satisfy a QM
			#			If a physician gets a "yes" for one patient, +0.30
			#			If a physician gets a "no" for another pt,   -0.70
				yes <- 1-avg
				no <- -1*avg
			
		#Output percentage difference from average of patients seen by each physician who satisfied QM
		print("Calculating quality measure ...")
		qm_phys <- t(sapply(phys, 
					  function (x){ 
					  				#Subset x to just eligible patients
					  				code <- x[1,"physcodeyear"]
					  				x <- x[x[,denominator],]
					  				#Recenter numerator logical vector around average
					  				x$numerator_avg <- ifelse(x[,numerator], yes, no)
					  					#Check if MD has any patients eligible for QM
					  					if(nrow(x[x[,denominator],])==0) return(as.numeric(c(code, 0 , NA)))
					  					if(nrow(x[x[,denominator],])>0){
					  						return(as.numeric(c(
					  							                  code, 
					  							                  nrow(x[x[,denominator],]), 
					  							                  sum(x$numerator_avg*x$patwt)/sum(x$patwt))))
					  					}		                  
					  		      }			
					 ))
		#Rename columns to be clear, convert physcodeyear to a character
		print("Outputting results ...")
		qm_phys <- data.frame(as.character(qm_phys[,1]), qm_phys[,2], round(qm_phys[,3], 3))
		names(qm_phys) <- c("physcodeyear", paste(numerator, ".pts", sep=""), numerator)
		return(qm_phys)
}

###########
#Section 3: Creating quality measures
###########
#Create physician level data.frames for manipulation and quality measure storage
	#Create NAMCS summary file on physician level
	namcs_phys <- namcs[!duplicated(namcs$physcodeyear),]
					
	################################################
	# 1) Blood pressure measurement in HTN patients#
	
	#Denominator: all visits from hypertensive patients
		#ICD-9 codes (from AMA-PCPI):
		#  401.(0|1|9), 40(2|3).(0|1|9)(0|1), 404.(0|1|9)[0-3], 
		namcs$htn.visits <- c(namcs$htn | diags(c("^1401(0|1|9)+", "140(2|3)(0|1|9)(0|1)", "1404(0|1|9)[0-3]")))
	#Numerator: all hypertensive patients who had their BP measured
		namcs$htn.measure <- c( !(namcs$bpsys %in% c(-9,999) | namcs$bpdias %in% c(-9,999)) & namcs$htn.visits)

	####Physician level calculation#####
	
	htnmeas_phys <- phys.level(numerator = "htn.measure", denominator = "htn.visits")
	namcs_phys <- merge(htnmeas_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")
	
	#################################################
	# 2) Blood pressure <140/90 in diabetes patients#
	
	#Denominator: all visits from DM patients by drug use or ICD-9 diagnosis (NCQA)
		#ICD-9 codes (from NCQA): 250.x, 357.2, 362.0, 366.41, 648.0 
		#Drug classes: 216 (alpha-glucosidase inhibitors), 314 (antidiabetic combinations), 215 (insulin), 282 meglitinides,
		#			   309 (misc agents), 213 (sufonylureas), 271 (thiazolidinediones)   
		#Exclusions: ICD-9: 256.4, 251.8, 962.0, 648.8 (PCOS, steroid induced DM, gestational DM)
		namcs$dm.visits <- c(
							(namcs$diabetes | 
							 diags(c("22050"), type="rfv") |
							 diags(c("^1250+", "^13572+", "13620+", "136641", "^16480+")) |
							 diags(c("216", "314", "215", "282", "309", "213", "271"), type="class")) &
							 !diags(c("^12564+", "^12518+", "^19620+", "^16488+"))
							)
	#Numerator: all hypertensive patients who had their BP measured
		namcs$htn.dm<- c( (namcs$bpsys < 140 & namcs$bpdias < 90) & namcs$dm.visits)

	####Physician level calculation#####
	
	htndm_phys <- phys.level(numerator = "htn.dm", denominator = "dm.visits")
	namcs_phys <- merge(htndm_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")

    ##################################################################
	# 3) Blood pressure <140/90 in ischemic vascular disease patients#
	
	#Denominator: all visits from IVD patients by ICD-9 diagnosis
		#ICD-9 codes (from NCQA): 411, 413, 414.(0|8|9), 429.2, 433-4, 440.1-2, 444, 445 
		namcs$ivd.visits <- c(
							 namcs$cebvd | namcs$ihd |
							 diags(c("^141(1|3|4)+", "^1414(0|8|9)+", "^14292+", "^143(3|4)+", "^1440(1|2)+", "^144(4|5)+")) 
							)
	#Numerator: all hypertensive patients who had their BP measured
		namcs$htn.ivd <- c( (namcs$bpsys < 140 & namcs$bpdias < 90) & namcs$ivd.visits)

	####Physician level calculation#####
	
	htnivd_phys <- phys.level(numerator = "htn.ivd", denominator = "ivd.visits")
	namcs_phys <- merge(htnivd_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")
	
	############################################
	# 4) Blood pressure <140/90 in all patients#
	
	#Denominator: all visits from adults patients, excluding DM and IVD
		namcs$all.visits <- !(namcs$ivd.visits | namcs$dm.visits)
	#Numerator: all patients with < 140/90 BP 
		namcs$htn.control.all <- c(namcs$bpsys < 140 & namcs$bpdias < 90)

	####Physician level calculation#####
	
	htnall_phys <- phys.level(numerator = "htn.control.all", denominator = "all.visits")
	namcs_phys <- merge(htnall_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")

	##############################
	#		Sum up the BP measures
	##############################	
	namcs_phys$bp.avg <- rowSums(namcs_phys[,c("htn.control.all", "htn.dm", "htn.ivd", "htn.measure")], na.rm=TRUE)/
						 rowSums(cbind(
						 		  as.logical(namcs_phys$htn.control.all.pts), as.logical(namcs_phys$htn.dm.pts), 
						 		  as.logical(namcs_phys$htn.ivd.pts), as.logical(namcs_phys$htn.measure.pts)
	  							))
	
	
	
	#*****CHF measures might not be appropriate****#
	############################################
	# 5) Beta blocker in CHF                   #
	
	#Denominator: (from AMA-PCPI)
		#ICD-9: 402.(0|1|9)1, 404.(0|1|9)(1|3), 428.0-1, 428.[2-4][0-3], 428.9
		#Exclusions: (only if didn't get beta-blocker) 427.8(1|9), 458.0-1, 458.2(1|9), 458.8-9, 493.00-92, 426.0,12,13  
	chf <- c("^1402(0|1|9)1", "^1404(0|1|9)(1|3)", "^1428(0|1)+", "^1428[2-4][0-3]", "^14289+")
	namcs$chf.visits <- c(namcs$chf | diags(chf))   
	 
	#Numerator: beta--blockers - 386, 274, 275
	namcs$chf.beta <- c(diags(c("274", "275", "386"), type="class") & namcs$chf.visits)		#Apply exclusions:
		beta.excl <- diags(c("14278(1|9)", "^1458(0|1)+", "14582(1|9)", "^1458(8|9)+", "^1493+", "^14260+", "14261(2|3)"))
		namcs$chf.beta.visits <- namcs$chf.visits & !(beta.excl & !namcs$chf.beta)
					
	chfbeta_phys <- phys.level(numerator = "chf.beta", denominator = "chf.beta.visits")
	namcs_phys <- merge(chfbeta_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")   							
	###########################################
	# 6) ACE/ARB in CHF                       #
	#Numerator: ACEI/ARBs: 042, 056
		namcs$chf.ace <- c(diags(c("042", "056"), type="class") & namcs$chf.visits)		#Apply exclusions: (from AMA-PCPI)
		ace.excl <- diags(c( "12776+", "^139(5|6)(0|2|8)", "^1403(0|1|9)1","^1404(0|1)(2|3)", "^14049(2|3)", 
							"^14251+", "^14401+", "^1584[5-9]+", "^1585(5|6)+", "^1586+", "174722", "^17885+", 
							"^20560+", "20568+", "^103995+", "^105498+", "^202[2-3]+"))  
		namcs$chf.ace.visits <- namcs$chf.visits & !(ace.excl & !namcs$chf.ace)
					
	chface_phys <- phys.level(numerator = "chf.ace", denominator = "chf.ace.visits")
	namcs_phys <- merge(chface_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")   							
   	############################################
	# 7) ACE/ARB in CAD with DM                #
	#Denominator: patients with CAD who also have DM or LVSD (from AMA-PCPI)
		namcs$cad.dm.visits <- c(
									(namcs$ihd |
									diags(c("^14140[0-7]", "^1414(8|9)+", "^1410[0-8]+", "14109[0-2]",  
									"^1412+",  "^1411[0-8]+", "^1413+", "20458(1|2)")) |
									diags("25150", type="rfv"))
									&
									(namcs$dm.visits)
						  		) 
	#Numerator: ACEI/ARBs: 042, 056
		namcs$cad.ace <- c(diags(c("042", "056"), type="class") & namcs$cad.dm.visits)
		namcs$cad.dm.visits <- namcs$cad.dm.visits & !(ace.excl & !namcs$cad.ace)
					
	cadace_phys <- phys.level(numerator = "cad.ace", denominator = "cad.dm.visits")
	namcs_phys <- merge(cadace_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear") 
	
	############################################
	# 7) Statins in CAD                        #
	#Denominator: patients with CAD (from AMA-PCPI)
		namcs$cad.statin.visits <- c(
									(namcs$ihd |
									diags(c("^14140[0-7]", "^1414(8|9)+", "^1410[0-8]+", "14109[0-2]",  
									"^1412+",  "^1411[0-8]+", "^1413+", "20458(1|2)")) |
									diags("25150", type="rfv"))
						  		  ) 
	#Numerator: Statins or other lipid-lowering agents
		namcs$cad.statin <- c(diags(c("173", "174", "241", "252", "316", "317"), type="class") & namcs$cad.statin.visits)
		statin.excl <- c(diags(c("^1995[0-1]+", "^19952(7|9)+")))
		namcs$cad.statin.visits <- namcs$cad.statin.visits & !(statin.excl & !namcs$cad.statin)
					
	cadstatin_phys <- phys.level(numerator = "cad.statin", denominator = "cad.statin.visits")
	namcs_phys <- merge(cadstatin_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")

	############################################
	# 7) Anti-thrombotic/aspirin in CAD        #
	#Denominator: patients with CAD (from AMA-PCPI)
		namcs$cad.asp.visits <- c(
									(namcs$ihd |
									diags(c("^14140[0-7]", "^1414(8|9)+", "^1410[0-8]+", "14109[0-2]",  
									"^1412+",  "^1411[0-8]+", "^1413+", "20458(1|2)")) |
									diags("25150", type="rfv"))
						  		  ) 
	#Numerator: Aspirin or other anti-thrombotics
		namcs$cad.asp <- c(diags(c("062", "211"), type="class") & namcs$cad.asp.visits)
		asp.excl <- c(diags(c("^1995[0-1]+", "^19952(7|9)+")))
		namcs$cad.asp.visits <- namcs$cad.asp.visits & !(asp.excl & !namcs$cad.asp)
					
	cadasp_phys <- phys.level(numerator = "cad.asp", denominator = "cad.asp.visits")
	namcs_phys <- merge(cadasp_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")

    ############################################
	# 8) Coumadin in AF                        #
	#Denominator: patients with AF with exclusions (from AMA-PCPI)
		namcs$af.visits <- c(
								diags(c("142731", "142732")) &
							   !diags(c("1394(0|2)+", "199602", "199671", "^20422+", "^20433+"))
							)
	
	#Numerator: Coumadin or related compound
		namcs$af.coum <- c(diags(c("262"), type="class") & namcs$af.visits)
				
	afcoum_phys <- phys.level(numerator = "af.coum", denominator = "af.visits")
	namcs_phys <- merge(afcoum_phys, namcs_phys, by.x="physcodeyear", by.y = "physcodeyear")
	
	#Average CVD measures
	namcs_phys$cvd.avg <- rowSums(namcs_phys[,c("chf.beta", "chf.ace", "cad.ace", "cad.statin", "cad.asp", "af.coum")], na.rm=TRUE)/						  rowSums(cbind(
						 		  as.logical(namcs_phys$chf.beta.pts), as.logical(namcs_phys$chf.ace.pts), 
						 		  as.logical(namcs_phys$cad.ace.pts), as.logical(namcs_phys$cad.statin.pts),
						 		  as.logical(namcs_phys$cad.asp.pts), as.logical(namcs_phys$af.coum.pts)
	  							))


######Finish up with imputation and save data

save(namcs, file="namcs_05-07.RData")
save(namcs_phys, file="namcs_phys_4-27-10.RData")

#Multiple imputation for missing data

#Drop raw data vars for analysis and imputation	
drop_vars <- c("age", "sex", "racer", "paytype", "totchron", "msa", "spec", "retypoff", "empstat", "owns", "emedrec", "edemog", "ecpoe", "ectoe", "eresult", "epnotes", "eremind", "prprvtr", "pccprod", "pccsat", "pccqoc", "pccpprof", "measpub") 
for (i in 1:length(drop_vars)) { namcs[,drop_vars[i]] <- NULL }

id_vars <- names(namcs)[c(1:7,13:14,16:17,19:59,78:98)]
nom_vars <- names(namcs)[c(60:69,71:75)]
ord_vars <- c("private", "numchron", "age.cat")

namcs_mi <- amelia(namcs, idvars = id_vars, nom = nom_vars, ord = ord_vars)
save(namcs_mi, file="namcs_mi_4-27-10.RData")

#Look at physician level data: quality measures by specialty

	#Create measure for all visits per physician
	namcs_phys$num_visits <- namcs_phys$htn.control.all.pts + namcs_phys$htn.dm.pts + namcs_phys$htn.ivd.pts
	namcs_phys$EHRscore <- namcs_phys$EHR + namcs_phys$demog + namcs_phys$escrip + namcs_phys$etest + namcs_phys$enotes + 
						   namcs_phys$result + namcs_phys$remind
	namcs_phys$anyEHR <- c(namcs_phys$EHR | namcs_phys$demog | namcs_phys$escrip | namcs_phys$etest | namcs_phys$enotes | 
						   namcs_phys$result | namcs_phys$remind)					   

scatter <- qplot(cvd.avg, bp.avg, data=namcs_phys, geom="point", color=as.factor(gp), size = num_visits, main="BP vs. CVD Quality Performance by Specialty", alpha=I(0.7), xlab="CVD Measure Performance", ylab="BP Measure Performance") + 
scale_colour_manual(name="Speciality", breaks=c(0,1), labels=c("Cardiologists", "Generalists"), values=c("red", "blue")) + scale_area(name="Number of Patient Visits", breaks=c(19, 31, 66), labels=c("<20", "20-30", "31+")) +
scale_x_continuous(breaks=c(seq(-0.6, 0.6, 0.2)), labels = c("-60%", "-40%", "-20%", "Average", "+20%", "+40%", "+60%")) +
scale_y_continuous(breaks=c(seq(-0.6, 0.2, 0.2)), labels = c("-60%", "-40%", "-20%", "Average", "+20%"), limits=c(-0.6, 0.25)) +
geom_hline(yintercept=0, linetype=2, legend=FALSE) + geom_vline(xintercept=0, linetype=2, legend=FALSE) +
theme_grey(base_size=14)


namcs_phys2 <- namcs_phys[!is.na(namcs_phys$EHR),]
scatter2 <- qplot(cvd.avg, bp.avg, data=namcs_phys2, geom="point", color=as.factor(EHR), size = num_visits, main="Blood Pressure  vs. CVD Medication Quality Performance by EHR Use", alpha=I(0.7), xlab="CVD Measure Performance", ylab="BP Measure Performance") + 
scale_colour_hue(name="EHR Use", breaks=c(0,1), labels=c("No EHR", "Any EHR Use")) + 
scale_area(name="Number of Patient Visits", breaks=c(19, 31, 66), labels=c("<20", "20-30", "31+")) +
scale_x_continuous(breaks=c(seq(-0.6, 0.6, 0.2)), labels = c("-60%", "-40%", "-20%", "Average", "+20%", "+40%", "+60%")) +
scale_y_continuous(breaks=c(seq(-0.6, 0.2, 0.2)), labels = c("-60%", "-40%", "-20%", "Average", "+20%"), limits=c(-0.6, 0.25)) +
geom_hline(yintercept=0, linetype=2, legend=FALSE) + geom_vline(xintercept=0, linetype=2, legend=FALSE) +
theme_grey(base_size=14)


arrange(scatter2, scatter, nrow=2, ncol=1)

